#include <iostream>
#include <fstream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TString.h"
#include "TEllipse.h"

using namespace std;

struct Step {
  int StepNum;
  float XPos,YPos,ZPos,Px,Py,Pz,edep;
};
struct PMT{
  int nhit;
  float x,y,z,px,py,pz;
  int isref;
};
struct Track{
  int pdgcode;
  float Edep;
  float PathLength;
};

void ana(){

  TFile* fread=new TFile("../outputFiles/cosmic1.root");
  if(!fread){std::cout<<"!!! No Input File !!!"<<std::endl;return;}

  TTree* tree=(TTree*)fread->Get("tree");

  int nEvent = tree->GetEntries();
  cout << nEvent <<endl;

  int EvNum;
  Step step;
  Track InitTrack;
  PMT JPMThit[100000],MPMThit[100000];
  int nPhoton,nJPMThits, nMPMThits;
  float EDEP;

  tree->SetBranchAddress("EvNum",&EvNum);
  tree->SetBranchAddress("Track",&InitTrack);
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nJPMThits",&nJPMThits);
  tree->SetBranchAddress("JPMThit",JPMThit);
  tree->SetBranchAddress("nMPMThits",&nMPMThits);
  tree->SetBranchAddress("MPMThit",MPMThit);
  tree->SetBranchAddress("EDEP",&EDEP);

  TH1F *hJPMTx=new TH1F("hJPMTx","hJPMTx",150,-150,150);
  TH1F *hJPMTy=new TH1F("hJPMTy","hJPMTy",150,-150,150);
  TH1F *hJPMTz=new TH1F("hJPMTz","hJPMTz",150,-150,150);
  TH1F *hMPMTx=new TH1F("hMPMTx","hMPMTx",150,-150,150);
  TH1F *hMPMTy=new TH1F("hMPMTy","hMPMTy",150,-150,150);
  TH1F *hMPMTz=new TH1F("hMPMTz","hMPMTz",150,-150,150);

  TH2F *hPMThit=new TH2F("hPMThit","hPMThit",150,-150,150,150,-150,150);
  TH2F *hPMThit_ref=new TH2F("hPMThit_ref","hPMThit_ref",150,-150,150,150,-150,150);

  for(int iev=0;iev<nEvent;iev++){

    tree->GetEntry(iev);

    std::cout<<"EvNum = "<< EvNum <<std::endl;
    std::cout<<"EDEP = " << EDEP <<std::endl;
    std::cout<<"track edep = " << InitTrack.Edep <<std::endl;
    std::cout<<"track PathLength = " << InitTrack.PathLength <<std::endl;

    std::cout<<"nJPMThits = " << nJPMThits <<std::endl;
    for(int iJhit=0;iJhit<nJPMThits;iJhit++){
      hJPMTx->Fill(JPMThit[iJhit].x);
      hJPMTy->Fill(JPMThit[iJhit].y);
      hJPMTz->Fill(JPMThit[iJhit].z);
      hPMThit->Fill(JPMThit[iJhit].x,JPMThit[iJhit].y);
      if( JPMThit[iJhit].isref==1 )  hPMThit_ref->Fill(JPMThit[iJhit].x,JPMThit[iJhit].y);
    }
    std::cout<<"nMPMThits = " << nMPMThits <<std::endl;
    for(int iMhit=0;iMhit<nMPMThits;iMhit++){
      hMPMTx->Fill(MPMThit[iMhit].x);
      hMPMTy->Fill(MPMThit[iMhit].y);
      hMPMTz->Fill(MPMThit[iMhit].z);
      hPMThit->Fill(MPMThit[iMhit].x,MPMThit[iMhit].y);
      if( MPMThit[iMhit].isref==1 )  hPMThit_ref->Fill(MPMThit[iMhit].x,MPMThit[iMhit].y);
    }
    
  }//for

  //-------------- Draw
  TCanvas *c1=new TCanvas("c1","c1",1200,900);
  c1->Divide(3,2);
  c1->cd(1);
  hJPMTx->Draw();
  c1->cd(2);
  hJPMTy->Draw();
  c1->cd(3);
  hJPMTz->Draw();
  c1->cd(4);
  hMPMTx->Draw();
  c1->cd(5);
  hMPMTy->Draw();
  c1->cd(6);
  hMPMTz->Draw();

  TCanvas *c2=new TCanvas("c2","c2",800,600);
  float pstep1=40.;
  float pstep2=pstep1*sqrt(2);
  TEllipse *eJPMT=new TEllipse(0.,0.,32.,32.);
  TEllipse *eMPMT=new TEllipse(-2*pstep1,0.,23.,23.);
  c2->Divide(2,1);
  c2->cd(1);
  hPMThit->Draw("colz");
  eJPMT->SetFillStyle(0);
  eJPMT->Draw();
  eMPMT->SetFillStyle(0);
  eMPMT->Draw();
  c2->cd(2);
  hPMThit_ref->Draw("colz");
  eJPMT->SetFillStyle(0);
  eJPMT->Draw();
  eMPMT->SetFillStyle(0);
  eMPMT->Draw();


}
